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1 Introduction 



As the predictions of QCD become increasingly precise and more high quahty data are available, 
the theoretical uncertainties associated with the analysis of the data are more and more often 
found to be dominant, and have thus come under increased scrutiny The determination 
of the QCD coupling as{Mz) from scaling violations of deep inelastic structure functions 
is perhaps the simplest and most fundamental example of this situation. In principle, deep 
inelastic scattering would be expected to provide a theoretically solid and experimentally clean 
way of determining as{Mz)- In practice, however, the situation is far from satisfactory |^, as 
revealed by the lack of stability of the value of as obtained through this procedure. 

The main source of difficulties can be traced to the fact that conventional extractions of 
from scaling violations of structure functions involve a simultaneous determination of parton 
distributions of the nucleon. The error on a^, therefore, gets tangled with the uncertainty on 
parton distributions, which is notoriously difficult to assess, and subject to a variety of sources of 
theoretical bias This is a consequence of the fact that the evolution equations of perturbative 
QCD 0], which govern scaling violations, can only be solved analytically in compact form by 
taking Mellin moments, which diagonalize the evolution operator. In practice, however, Mellin 
moments are not directly measurable, since their determination would require the availability of 
data taken at arbitrarily large center-of-mass energy. The problem is circumvented by solving 
directly the momentum-space Altarelli-Parisi equations, which, however, usually entails the 
construction of a parton parametrization. 

Such an undertaking thus runs into two difficulties. First, a parton parametrization is 
routinely constructed by assuming a functional form for parton distributions, and fitting the 
corresponding free parameters. The choice of a functional form, however, introduces a theoret- 
ical bias, which is potentially large, and whose size is very difficult to assess For example, 
choices of functional form based on "QCD expectations" may well give a poor representa- 
tion of unexpected or unconventional behaviors of the distributions. Furthermore, any given 
parametrization constrains the form of the distribution near and beyond the boundary of the 
region where data are available: whenever the data are not very precise, it can be seen 
that the results obtained for observable quantities may depend significantly on the form of the 
parametrization. 

Second, it is quite difficult to estimate precisely the error on parton distributions, and 
even harder to propagate this error to quantities which depend upon them. Indeed, a generic 
observable (such as a Mellin moment, or a cross section) is a functional of parton distributions, 
which in general depends on it in a complicated, nonlocal way. Thus, on the one hand it is 
difficult to extract the errors on parton distributions from errors on the data, on the other 
hand it is difficult to propagate errors on parton distributions to quantities calculated from 
them. At the very least, the full information on experimental errors and correlations should 
be used, but it may also turn out that conventional techniques of linear error propagation fail. 
These difficulties could only be overcome if one were able to obtain a reliable and unbiased 
representation of the probability measure in the functional space of parton distributions , as 
determined by the available experimental information. 

In this paper, we approach the determination of as in a way which bypasses these difficulties, 
by combining two techniques which have recently been proposed and implemented in the context 
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of the analysis of DIS data: the use of evolution equations for truncated Mellin moments of 
parton distributions P-p^ , and the determination of the probability measure in the space of 
deep inelastic structure functions by means of neural networks [|l^ . 

Truncated moments are defined as Mellin moments with the integration range restricted to 
the subset xq < x < 1 of the kinematically allowed range < x < 1. Unlike ordinary Mellin 
moments, they are measurable quantities, since the small x region is excluded. The evolution 
equations which they satisfy turn out to be a system of coupled ordinary linear differential 
equations, which can be truncated to finite accuracy to yield a closed system. The evolution 
equations for truncated moments have been derived, and explicitly solved at next-to-leading 
order in the nonsinglet |^ and singlet [0 sectors. The numerical accuracy of the approximate 



solutions has been studied and improved upon by means of suitable techniques [|TT|, and is 
equal or better than that of standard numerical solutions of Altarelli-Parisi equations [|l|. 

Given an experimental determination of truncated moments at more than one scale, one 
can then determine as directly by comparing computed and measured scaling violations. In the 
nonsinglet sector, there is only one independent structure function, say F2, and scaling violations 
are fully determined by and an initial condition, given by the values of the truncated moments 
of F2 at a reference scale: therefore one does not have to use a parton parametrization. 

In practice, however, a sufficiently accurate determination of truncated moments of the 
nonsinglet F2 and associated errors and correlations cannot be obtained by simply binning 
and summing data points: the coverage and precision of the data are not sufficient for this 
purpose. Rather, in order to fully exploit the available experimental information, it is necessary 
to construct a parametrization of the structure function in the kinematic region where data 
are available. This clearly leads back to the same difficulties encountered when constructing 
a parton parametrization. The problem, however, can now be overcome, since an unbiased 
determination of the probability density in the space of structure functions, based on available 
experimental data, has been constructed in Ref. fT^ . 

The representation of the probability density given in Ref. |JT^ takes the form of a set of 
neural networks, trained on an ensemble of Monte Carlo replicas of the experimental data, 
which reproduce their probability distribution. The parametrization is unbiased, since neural 
networks do not require the choice of a specific functional form, and it interpolates between data 
points, imposing smoothness constraints in a controllable way. Information on experimental 
errors and correlations is incorporated in the Monte Carlo sample, so that errors on truncated 
moments and correlations between them can be determined without having to use linearized 
error propagation. The parametrization combines all the available experimental information, 
so that statistical errors are minimized, but it also correctly reproduces the loss of accuracy 
incurred when extrapolating outside the kinematical region where data are available. 

Hence, we can obtain an unbiased evaluation of truncated moments, and then use this to 
determine directly without any further assumptions. We end up with a determination of as 
where all sources of uncertainty related to the method of analysis are under control, and the 
only theoretical uncertainty is related to the lack of knowledge of the anomalous dimensions 
beyond next-to-leading order. This gives us a bias- free determination of a^, and simultaneously 
illustrates the power of a method of analysis based on the direct knowledge of a probability 
density in a space of functions. 

This paper is organized as follows: in section ^ we will review the method of truncated 
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moments emphasizing the issues of numerical accuracy which are relevant to the rest of 



the paper; in section |^ we will introduce the neural network parametrization of DIS structure 
functions developed in Ref. [|12|; section ^ contains the details of our fitting procedure and 



our result for the strong coupling as{Mz): we explain our data selection, our choice of fitting 
architecture, our error estimate and the consistency tests that we performed; finally, section § 
summarizes our conclusions and discusses possible extensions of our work. 

2 Truncated moments of parton distributions 

We determine the strong coupling from the scale dependence of the nonsinglet deep inelastic 
structure function 

F,'''ix,Q') = F^ix,Q') - F,\x,Q') . (2.1) 
In the DIS scheme -^2^'^ is expressed directly in terms of the nonsinglet combination of 



quark distribution, according to 

ef |(7i(x,Q^) + q^ix,Q'^) 

J p—n 



Ff^(x,Q2) = Y.e^ [q.{x,Q') + q,{x,Q')] _^ . (2.2) 



i=l 



The scale dependence of F2 is thus given by the evolution equation for nonsinglet quark distri- 
butions, henceforth denoted simply by g(a;,/i^). 



Jx y \y J 



/^^ A ?(^. f^') = ^ -P ( «.(/^') 1 liy^ f^') . (2.3) 



where P{z, as) is the appropriate evolution kernel. 

2.1 Evolution equations for truncated moments 

Consider now the truncated Mellin moment of the quark distribution 



qn{xQ,^^)= [ dx x"- ^g(x,/i^) . (2.4) 



The evolution equation for truncated moments is easily found to be 



/^^^ qnixo, n') = t dy f'-'q{y, ^^) G„ ( ^; a,(/i^) ) , (2.5) 

where Gn is the truncated moment of the splitting function 

Gn{x,as) = f' dzz''-'P{z,as) . (2.6) 

Note that the function Gn{x, as) is analytic on the real axis in the interval < x < 1, while it 
has integrable singularities at x = 1 as a consequence of the presence of + distributions in the 
kernel P{z, as) at all perturbative orders. 
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We can turn eq. (|2.5| ) into a set of coupled linear differential equations by expanding 
G{xo/y, as) in powers of y around y = 1 as 



Gn f-, = E ^ g'^\xo. as) {y - 1)^ . (2.7) 

Because the singularities of Gnixo/y, as) aX y —* xq are integrable, one can substitute eq. ( |2.7|) 
into eq. (p.5|) and integrate term by term, obtaining a convergent series. It is then legitimate 
to truncate the sum in eq. ( |2.7| ) at a finite order p = M, with the result 



Gn[j,as^=f:4f{xo,as)y' + 0{y-l)''+' , (2.^ 

where 



i':^\xo,as) = t ^~Kff^7'''^ - (2-9) 
k=p P\k~p)\ 



The evolution equation for truncated moments then becomes 

/^^;A?n(a^o,/i^) = ^ E 4*^)(xo,as) gn+p(a;o,/i^) , (2.10) 



' 27r 



p=0 



and is thus given by a set of coupled linear differential equations. The anomalous dimension 
matrix c^^^(xo, Og) couples each truncated moment g„ to all truncated moments qk with k > n. 
The evolution equation for the k-th moment can be made arbitrarily accurate by taking M to 
be sufficiently large. 

2.2 Solution of the evolution equations 

The coupled evolution of the moments qk, with n < k < n + M, can be solved to determine the 
value of qn at different scales, if the system in eq. ( |2.1U| ) closes, i.e. if it is legitimate to include 



a decreasing number of terms in the evolution equation of increasingly higher moments. In 
Ref. [§] it was shown that this is indeed the case, because the matrix elements of c^^^(a;o, ctg) 
decrease rapidly as one moves away from the diagonal, so that the accuracy of the truncated 
evolution kernel actually improves with increasing order of the moment, even if one less term 
is included when the order is increased by one unit. 

Having chosen the values for n and M, we can introduce the simplified notation 

Cki{as) = 4ti~k^''\xo, as) {n + M > I > k > n) , ,^ 
Ckiias) = {n<l<k<n + M), ^' ^ 

in terms of which the system of equations to be solved reads simply 

f^^:r^1k{xo,fi^) = — E Cki{as) qi{xo,fx'^) , (2.12) 

d/x2 27r ^ 
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with n < k < n + M. Expanding Cki{ois) in powers of as to next-to-leading order, one can 
now solve eq. ( |2.12| ) by standard perturbative methods. The task is considerably simplified by 
the fact that the matrix of coefficients, Cki{as), is upper triangular, so that the eigenvalues 
are exactly the diagonal elements Ckk^ and the eigenvectors are easily computed by means of a 
recursion relation. The explicit solution is given in Ref. 0. 

In practice, the numerical implementation of the solution may be problematic if a high 
accuracy is required for low moments. Indeed, it turns out that when xq = 0.1 the value of 
M which is needed to achieve an accuracy at the few percent level on the evolution of the 
ra-th truncated moment is typically of order M > 100 for the lowest moments, n = 1, 2 (when 
n < 1 the convolution integral in eq. (p.5|) diverges) , even though it rapidly decreases for higher 
moments. Such a large value of M leads to numerical difficulties, due to the fact that the 
matrix elements Cki{as) become very small when I » k. 



2.3 Improved solution 

An improved solution to the evolution equation, which overcomes the problems related to 
the slow convergence of the expansion of the evolution kernel for low moments, can however be 



derived [jTT| . The origin of the problem is easily tracked to the logarithmic singularities left over 
in the function Gn{xo/y,ag), for y — > xq, as a consequence of the presence of -|- distributions 
in the evolution kernel P{z,as), for z ^ 1. For the lowest moments, the integration over y 
is dominated by small values of ?/, y ~ Xq, where the function G„ is in turn dominated by 
the logarithmic singularity, and not well approximated by a low-order truncation of its Taylor 
expansion around y = 1. 

To solve the problem, one can integrate by parts the right hand side of eq. ( ^.51 ), obtaining 

/ dy y'^~^ q{y,fi'^) (— ,as) = \Gn{xo,y;as) y""^ g(z/, 

' dy G^ixo,y;as)-^ (y^~'qiy,Q')) , (2.13) 
'XO ay ' 



1 

XO 



where 



Gn{xo,y;as) = [ dzGn(—,as) . (2.14) 

J XO \ Z 



In the definition of Gn the freedom to fix the lower limit of integration, which is independent 
of y, has been exploited, so that the function G„ satisfies G„(xo, Xq; as) = 0, as well as 

Gn{xo,y;as) = y Gni^,asj - xo Gn-ii^,as] ; (2.15) 



furthermore, the coefficients of the Taylor expansion of G„ around y = 1, defined as in eq. 
satisfy, in obvious notations, 

^^")(xo, as) = g'p'Mxo, as) ■ (2.16) 

The integration region ?/ ~ xq is now suppressed, and a faster rate of convergence when G„ 
is expanded around y = 1 is expected. Repeating the procedure that leads to eq. ( p.l2|) one 
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gets 



11 



27r 



p=0 P- 



M 



1=1 



(2.17) 



The coefficients Cni{ots) are given in eq. ( p.ll|) . The ffist term on the right hand side of 
eq. ( |2.17| ) is responsible for the faster convergence of the expansion to the exact result: it 
vanishes as M — > cxd, as it must, since it reconstructs G„(a;o,a;o; ds) = 0. 

Eq. (|2.17|) , however, cannot be used directly: as the right-hand side of the equation depends 
on the value of g(xo,/i^), we do not get a closed system of equations. This obstacle can be 
circumvented by noting |jll| that q{xo, n'^) itself can be expanded over a basis of truncated 
moments, and, specifically, a finite set of them is sufficient to give an accurate determination of 
q{xo, fi"^). To see this, expand the quark distribution in Taylor series about a selected point, say 
Uo = {1 + Xq)/2 (note that ?/o > (1 + Xq)/2 is not allowed since y = 1 is an essential singularity 
of all parton distributions) 



g(x,/i^) -- 

Truncated moments are then given by 



k=l 



(2.18) 



k=l 



where 



f3jk = dx x^ '^{x - yo) 



Xo 



(2.19) 



(2.20) 



The infinite matrix f3jk can be approximated by truncating it to a square N x N matrix (3jk, 
which is easily computed and inverted. One finds then an approximate expression for q{xQ, fi"^), 
given by 

N N ^ 

g(xo,/x') = EE/?fc"/?.(^o,/i')(xo-Z/o)'-' • (2.21) 
fc=ii=i 

The error introduced by the truncation of the matrix jSjk has been studied, as a function of A^, 



in Ref. |]TT|. It turns out that a satisfactory accuracy (at the 0.1% level) can be reached already 
for ~ 5, while in practice the inversion of the matrix jSjk becomes numerically difficult only 
for > 10 (as easily verified, the matrix is ill-conditioned). The method is thus viable and 
does not lead to loss of accuracy. 

Combining eq. ( p.l7|) and eq. ( p. 21 ), we get the final form for the evolution equation, which 
replaces eq. ( p.l2|) , namely 



fx 



d 
diJ? 



qk{xQ,^J? 



2tt 



n+M 



^ [Cui{as) + D 



kl ' 



qiixo^ij^) 



(2.22) 



l=n 
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where 



k-l 



M 



J2- 9pixo,as)ixo-iy 



N 



Once again, eq. ( p. 221 ) can be solved by expanding Cki{as) and D]^i'{as) to next-to- leading 



Lr=l 



,r-l 



(2.23) 



order in ag, and then using well-known perturbative methods. A compact expression for the 



solution is given in Refs. |TD|,|TT|. The price to pay for the faster convergence of eq. ( [^.22D is 
that the evolution matrices are no longer triangular. This implies that all truncated moments 
with n > 1 must be used as initial conditions in order to compute the evolution of any moment. 
As a consequence, diagonalization has to be performed numerically rather then analytically. 
In practice, however, this does not lead to any numerical difficulties: as shown in Ref. an 
accuracy at the few percent level on the evolution can be achieved for any truncated moment 
with M ~ 10, a size at which numerical diagonalization is not difficult. It is important to 
notice that all the approximations that have been introduced in order to make the problem 
amenable to a numerical solution are under theoretical control, and the accuracy achieved can 
be precisely estimated and can be improved upon if the need arises. 



3 Neural network parametrization of deep inelastic struc- 
ture functions 

Ideally, a parametrization of structure functions must incorporate all the information contained 
in the experimental measurements, i.e. their central values, their statistical and systematic er- 
rors, and their correlations; furthermore, it must interpolate between them without introducing 
any bias. Several approaches to this problem have been proposed, in the context of fitting par- 
ton distributions: ^^-minimization of a fixed functional form coupled to error propagation 
through the covariance matrix, and various improvements thereof [|1^]; projection over bases of 
orthogonal polynomials |15]; or Monte Carlo sampling coupled to Bayesian inference Here 
we will follow the method of Ref. , where an unbiased extraction of the probability measure 
in the space of structure functions is performed, based on a coordinated use of Monte Carlo 
generation of data and neural network fits. 



3.1 Experimental data 

Since we are interested in the nonsinglet structure function F2, we need a simultaneous mea- 
surement of this structure function for proton and deuterium targets. We will use the data of 
the New Muon Collaboration (NMC) ^ and of the BCDMS (Bologna-CERN-Dubna-Munich- 
Saclay) Collaboration |l7|, which provide a simultaneous determination of the proton and 
deuteron structure functions in the same kinematic region, and provide the full set of cor- 
related experimental systematics for these measurements. Earlier data from SLAC are not 
competitive with these in terms of accuracy and kinematic coverage. The more recent HERA 
data are available in a much wider kinematic region, but only for proton targets. Another set of 
joint proton and deuteron measurements was performed by the E665 Collaboration |jl8|. These 



data however are mostly concentrated at low Q (and low x), and thus are not relevant for 
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Figure 1: NMC and BCDMS kinematic range. 



perturbative QCD applications. The kinematic coverage of the data which we include in our 
analysis is displayed in Fig. |T]. Altogether our parametrization is based on N^at = 552 exper- 
imental points for the nonsinglet structure function, obtained as difference of pairs of proton 
and deuteron data according to eq. Henceforth, F^'^^^'^ will denote the i-th data point 



3.2 Probability measure in the space of structure functions 

The experimental data give us a probability measure in an A^^rfaf-dimensional space, assumed 
to be multigaussian. In order to extract from it a parametrization of the desired structure 
function we must turn this measure into a measure ^[-^2] in a space of functions. Once such a 
measure is constructed, the expectation value of any observable J-' [-^2(0;, Q^)] can be found by 
computing the weighted average 



F2{x,Q^) V[F, 



(3.24) 



Errors and their correlations can also be obtained from this measure, by considering higher 
moments of the same observable with respect to the probability distribution. 

The determination of an infinite-dimensional measure from a finite set of data points is an 
ill-posed problem, unless one introduces further assumptions. For instance, one may assume 
a fixed functional form [|14|, in which case the problem is reduced to the determination of a 
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finite set of parameters, or one may project over a basis of orthogonal polynomials [1^, and 
choose a particular truncation of the a priori infinite set of basis functions. In the approach 
of Ref. |T^, neural networks are used as interpolating functions, so that the only assumption 
is the smoothness of the structure function. Neural networks can fit any continuous function 
through a suitable training; smoother functions require a shorter training and less complex 
networks. Hence, an ideal degree of smoothness can be established on the basis of a purely 
statistical criterion (such as having reached a satisfactory goodness-of-fit) without need for 
further assumptions. 



3.3 Fitting strategy 

The construction of the probability measure is done in two steps: first, a set of Monte Carlo 
replicas of the original data set is generated. This gives a representation of the probability 
density V[F2] at points {xi,Qf) where data are available. Then, a neural network is fitted to 
each replica. The ensemble of neural networks gives a representation of the probability density 
for all {x,Q'^): when interpolating between data the uncertainty will be kept under control 
by the smoothness constraint, but it will become increasingly more sizable when extrapolating 
away from the data region. 

The k = 1, . . . , Nrep replicas of the data are generated as 



rp{art)(k) _ /i , Jk) ^ 
- U + Tjjv ai^N) 



rp{exp) Z^g t,aJt,a p(exp) (A:) 

i "T j^gg "T ' l,S i,S 



(3.25) 



where f/'"^^^ = F^ix^Ql) are the original data, fi^a are the experimental systematic errors, given 
in percentage, ai^N is the experimental normalization error, ai^s is the experimental statistical 
error, while rj^^-*, r^^} and are univariate gaussian numbers. When two data points share 
a correlated error or normalization, the same gaussian number is used. In Ref. |T^, a set of 



Nrep = 1000 replicas of this form has been generated, and it has been verified that central 
values, errors and correlations of the original experimental data are well reproduced by taking 
the relevant averages over a sample of this size. The kinematic bound ^2(1,(5^) = has been 
implemented by adding a number of artificial data points at x = 1 with carefully adjusted 
errors. 

Each set of generated data is fitted by an individual neural network. A neural network |T^ 



is a function of a number of parameters, which fix the strength of the coupling between neurons 
and the thresholds of activation for each neuron. The architecture of the network is chosen 
to be redundant for the given problem, i.e. the number of parameters is rather larger than 
the minimum needed to get a good fit. The parameters are then fitted by minimizing an error 
function. The fit is done using the technique of learning by back-propagation, whereby the data 
are repeatedly shown to the network until satisfactory learning is achieved. The error function 
is 



E(^) = y ^ . ^, (3.26) 

where jg ^j^e prediction for the i-th data point from the net trained on the k-th replica 

of the data. Use of this error functions, which only includes statistical errors, ensures [1^ 
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that statistical errors on adjacent data points are correctly combined, whereas the correlated 
systematics are reproduced when averaging over the Monte Carlo sample. 



3.4 Results and validation 



In Ref. |T^, an independent set of neural networks was trained on the non-singlet combination 



f/"^-*. This procedure is advisable because the structure functions are roughly of the same 
size, while their difference is typically by a factor 10 smaller. The length of training is fixed 
by studying the behaviour of the error function E^^\ defined as in eq. ( p.26| ) for the neural net 
fitted to the central experimental values, and asking that E^^yNdat stabilizes to a value close 
to one. 

A number of checks is then performed in order to make sure that an unbiased representa- 
tion of the probability density has been obtained. For instance, it has been verified that the 
covariance of two data points computed from the Monte Carlo sample of nets is on average very 
close to the corresponding covariance matrix element of the data. Since correlations of the data 
are entirely due to systematics, this indicates that the systematics is correctly reproduced. 

It turns out that the average standard deviation for each data point computed from the 
Monte Carlo sample of nets is substantially smaller than the experimental error. This could 
be due to the fact that the network is combining the information from different data points by 
capturing an underlying law, or that it is introducing a smoothing bias. In Ref. [|1^] a statistical 
indicator which distinguishes the two cases was constructed: define 



{net)(k) 



i=l 



a. 



(exp)^ 



(3.27) 



i. e. a modified error function where the prediction of the k-ih net for the i-ih. point is compared 
to the central experimental value rather than to the /c-th replica. The desired indicator is 



n 



rep 



(3.28) 



rep 



where {)rep denotes the average over the ensemble of replicas, and E, E are respectively defined 
by eqs. (|3.26| ) and ( p.27|) . One can then show that, in the limit in which the error computed from 
the Monte Carlo sample is much smaller than the experimental one, 7?. ~ 1/2 if an underlying 
law is captured by the net, but 7^ > 1 if the error reduction is due to a smoothing bias. The 
Monte Carlo sample of nets of Ref. |]T2[ has TZ = 0.58. 
The final set of neural networks F, 



Jnet){k) 



provides a representation of the probability measure 
in the space of structure functions, which can be used to estimate any functional average, defined 
as in eq. ( ^.24| ), using 



Nr. 



J2 T F2("^*)('=)(a;,Q') 



rep 



(3.29) 



In particular, the average, variance, and covariance of truncated moments computed using 
the Monte Carlo sample will give us a determination of values, errors and correlations of the 
measured truncated moments. 
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Table 1: 

Errors and correlations for various truncated moments with xq = 0.03 and = 20 GeV^. 
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Table 2: 

Errors and correlations for the fourth truncated moment with xq = 0.03 at various scales 
20 < g2 < 70 GeV2. 

4 Determination of 

Using the neural parametrization of structure functions, we can compute directly experimen- 
tal values of truncated moments at any scale in the region of the data. Because the neural 
parametrization retains all the experimental information, and specifically it allows a determi- 
nation of errors and correlations, we can view such values as direct experimental determinations 
of the truncated moments. The value of ag can then be extracted from scaling violations of 
truncated moments. Specifically, given the full set of truncated moments at a reference scale, 
the value of any truncated moment at any other scale is predicted in terms of as- Hence, we 
can obtain a determination of as by comparing this prediction with the data, for any choice 
of moment and scale. Even though such predictions are clearly correlated, it is to be expected 
that the use of a larger number of moments or scales will in general lead to a more precise deter- 
mination of ag] the accuracy should, however, saturate for a large enough number of moments 
and scales. Errors and correlations for typical truncated moments are displayed in Tables ^ 
and 0. Correlations are quite large, so it is conceivable that an optimal fit may be obtained 
with a relatively small number of moments. 

It is clear that to achieve a best fit for the strong coupling we must optimize the choice of 
the fitting procedure and parameters. We now turn to a detailed explanation of our choices for 
the truncation point of Mellin moments, for the range and number of scales Q^, for the number 
of moments to be used in the approximate evolution equation, the number of moments that 
should be treated as free parameters, and the effect of variations of these choices. We will then 
give our best fit with its associated statistical error. Finally, we will address the known sources 
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of theoretical error. 



4.1 Choice of truncation point and fitting range 

The choice of values of {xq, Q"^) and of the moments to be used for the extraction of is 
determined by the kinematic coverage of the data (see Fig. p, as reflected by the errors on 
the moments. In particular, at low there are no large x data, while at high there are 
no small x data. In view of the fact that the use of truncated moments allows one to exclude 
the small x region, and that at low power corrections can be sizable, it is convenient to 
consider only the large region. Also, the large x extrapolation is severely constrained by 
the kinematic bound -^2(1, Q^) = 0, whereas the small x behaviour is essentially unknown pO| . 
A reasonable cut which ensures a good coverage of the large x region is > 20 GeV^. By 
choosing a high enough truncation point, xq > 0.3, it would be possible to compute accurately 
truncated moments up to the highest available scale ~ 200 GeV^. Such a choice, however, 
is not advisable, because the bulk of the data would then be excluded. In fact, as Xq is raised, 
correlations between truncated moments rapidly increase. The value of xq should thus be chosen 
as small as possible, in order not to loose information. Requiring xq to be as low as xo = 0.03 
imposes then a cut < 70 GeV^. Smaller values of xq would require extrapolation. 

We take thus as a basehne choice xq = 0.03, 20 < < 70 GeV^. With this choice, mo- 
ments 1 < n < 8 have errors below 10%, and correlation coefficients are below 0.9, as long 
as neighbouring moments are avoided, and one does not consider more than three (logarith- 
mically) equally spaced scales in the available range. Correlations between moments are 
not significantly reduced by further lowering xq, which would anyway require extrapolation. If, 
however, xq is raised to Xq = 0.1, then correlations between neighbouring moments are greater 
than 0.98, and only correlations between moments differing by more than two orders are below 
0.9. Because correlations between the same moment evaluated at different scales cannot be 
reduced without enlarging the range, no more than three scales should be used if we wish 
to keep such correlations below 0.9. We will later consider variations of xq, the range and 
the number of scales about this choice. 



4.2 Evolution equation 



As discussed in section ^ the scale dependence of any truncated moment qn{xo,Q'^) can be 
determined to any required accuracy from the knowledge of a finite set of truncated moments 
at a reference scale Qq. The result has the form 



M 



(ln{xo,Ql)= Mnp{xo;Ql,Ql,as) qp{xo,Q, 



0) ; 



(4.30) 



p=n„ 



where the evolution matrix Mnp, explicitly given in Refs. |T^,|rT|, is determined as a perturbative 
expansion in ag. 

Given a measurement of truncated moments q^^^ixQ, Q"^) at more than one scale, the value 
of as can be determined by minimizing 

n,i m,j 



X 



(4.31) 



13 



n 


Xq = U.i 


^ n no 
Xq = U.Uo 


Xq = U.Ui 


o 

z 


0.091 ± 0.047 


0.085 ± 0.070 


0.089 ± 0.080 


6 


0.100 ± 0.024 


O.lOo ± 0.030 


0.106 ± 0.031 


A 

4 


0.113 ± 0.019 


0.115 ± 0.019 


0.115 ± 0.019 


o 


u.izz Hz u.uio 


n 1 9Q -u n ni 
u.izo Hz u.uio 


fi 1 9Q -u n ni 

U.lZO ZL U.UIO 


6 


0.127 ± 0.014 


0.127 ± 0.014 


0.127 ± 0.014 


7 


0.129 ± 0.015 


0.129 ± 0.014 


0.129 ± 0.015 


8 


0.129 ± 0.016 


0.129 ± 0.016 


0.129 ± 0.016 


9 


0.129 ± 0.018 


0.129 ± 0.018 


0.129 ± 0.018 



Table 3: 

Fits of as{Mz) from the evolution of a single moment. 



where Vni-mj is the covariance matrix for moments q^^{xQiQl), q^^{xQ,Q'j). Of course, the 
result for ag should be independent of the choice of initial scale Qq. 

If knowledge of moments qn{xo, Qq) with n^m < n < M is needed in order to obtain the 
desired accuracy in the solution of the evolution equation, then in principle all these moments 
should be treated as free parameters in the minimization, along with the value of as, so the sum 
over n, m in eq. (|4.31|) should run from rimin to M. In practice, however, for any reasonable value 
of Xq (say, Xq ~ 0.1), the evolution of truncated moment qn{xQ, Qq) is "almost diagonal", in the 
sense that it receives only a small correction from moments qm{xo,Qo)> with m ^ n. Because 
of this, and because of the large correlations between different moments, we may perform the 
minimization while only including a subset of moments in the sum over n, m in eq. ( |4.31|) , and 
treating only such moments as free parameters. This issue is discussed in detail in the next 
subsection. 



Throughout this section, evolution will be performed using the improved method of Ref. |Tl 
discussed in section 2.3. Specifically, we will take rimin = 1, with twelve moments included in 
the evolution equation (M = 11), while the reconstruction of q{xQ,Ql) according to eq. (|2.21|) 
will be performed with = 6. This ensures an accuracy of order 0.1% on the evolution 
of the second {n = 2) truncated moment, rapidly improving as n increases. For as, we use 
the solution of the next-to- leading order renormalization group equation to express ^^(Q^) in 
terms of as{Mz), which is then directly taken as a free parameter. The number of active flavors 
varies by one unit at each quark threshold, and the coefficients of the (3 function are matched 
according to the Marciano |21|| prescription. Dependence on quark thresholds will be studied 
as a source of theoretical uncertainty in section 



4.3 Choice of fitted moments 

Having fixed the number of truncated moments to be included in the evolution equation, one still 
has to choose which ones should be treated as free parameters in the minimization procedure, 
and which ones should be fixed at the experimental central value. The simplest possibility is 
to fit only one moment at a time. The results for as obtained in this case are displayed in 
Table ||, with the default choice of scales, and three different choices of truncation point xq. 
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0.03 


Fitted moments 




3+4 






0.137 ± 0.011 


2+4 






0.140 ± 0.010 


3+5 






0.136 ± 0.011 


4+6 






0.131 ± 0.012 


5+7 






0.128 ± 0.012 



Table 4: 

Fits of as{Mz) from the evolution of a pair of moments. 



The table shows that the uncertainty in the determination of function of the order of 

the moment n, has a minimum around n ^ Q. The presence of a minimum can be understood 
as a consequence of the fact that the most accurate data have large x ~ 0.5, but there are no 
data with x > 0.75; furthermore, the anomalous dimension vanishes at low n ~ 1 (it vanishes 
exactly at n = 1 when xq = 0), and increases monotonically in modulus as n increases. Hence, 
low moments lead to a less precise determination both because they are dominated by small 
X data and because their scaling violations are weaker; for high enough n, on the other hand, 
there is loss of precision due to the extrapolation in the very large x region. In particular 
(see also below) moments with n > 9 probe a region of x where elastic contributions to the 
cross-section start being relevant, and become rapidly unreliable as n increases. 

It is interesting to observe that all the determinations of Table ^ are compatible within 
errors. Also, it is interesting to observe that, for the values of n which give a reasonably precise 
determination of as, there is little or no dependence on the choice of Xq of both the central value 
and the error. The error on any of these individual determinations is however much larger than 
the error obtained from existing fits to these data |^,^: a more accurate determination can 
only be obtained by combining the information from different truncated moments, i.e. fitting 
more than one truncated moment at a time. 



When including more than one moment in the sum of eq. ( [4.3 1|) , the issue of correlations 



between moments becomes important. The impact of correlations is illustrated by the results 
obtained fitting a pair of moments, displayed in Table |. Because of large correlations, in each 
case the value of as turns out to be larger than either of the values obtained from each of the 
two moments individually. 

The presence of large correlations entails two distinct problems when performing a fit where 
several moments are simultaneously fitted. The first problem is that, as the elements of the 
covariance matrix Vij —>■ 1, the matrix becomes singular, in that all eigenvalues but one vanish 
in the limit. Hence, when correlations are large, several eigenvalues become very small and the 
inversion of the covariance matrix is numerically problematic. The second problem is that, even 
if the covariance matrix is accurately inverted, when correlations are large, off-diagonal terms in 
the GQ- ( |4.31|) may dominate over diagonal ones, thereby leading to generally unreliable and 
often pathological results. In particular, it may turn out that the best-fit values of all moments 
differ from the measured values by several standard deviations. An example of this pathological 
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Figure 2: 

Pathological best fit of the third moment in the presence of off-diagonal instabilities. 

situation is displayed in Fig. which shows the best-fit behaviour of the third moment {n = 3) 
in a fit of moments 2 + 3 + 4 + 5 + 6, chosen as a representative case in which moments are over- 
correlated because too many neighbouring moments are fitted.^ Even though this situation is 
in principle acceptable, in practice such results are unreliable because the minimum of the 
is obtained by a fine-tuned balance of diagonal and off-diagonal contributions, which cannot be 



trusted whenever errors and correlations are known with limited accuracy p5[. We will refer 
to this situation as off-diagonal instability of the fits. 

Because of these problems, there is a trade-off in accuracy when including more than one 
moment in the fit: the addition of more moments brings in new information, thus reducing 
statistical errors, but for a large enough number of moments it is impossible to accurately invert 
the covariance matrix, and fits are spoiled by instabilities. Whereas it is possible to keep the 
covariance matrix inversion under control by improving the numerical accuracy of calculations, 
the off-diagonal instability depends on the accuracy in the experimental determination of errors 
and correlations, and it is unavoidable. We have thus made sure that the covariance matrix is 
inverted to an accuracy which is by several orders of magnitude greater than the knowledge of 
its individual matrix elements. Then, we have tested for off-diagonal instabilities by flagging 
all results in which, at the best flt, one or more moments differ by more than one standard 
deviation from the experimental central value at any of the fltted scale, and discarding flts for 



^Note that even though this pathological situation is reminiscent of the effect which is found when linea rizing 
correlated normalization errors it is a distinct problem. To check this, we have replaced in the eq. ( 4.31 ) 



all Qn with In g„ , and computed the covariance matrix accordingly. In such case normalization uncertainties 
decouple, yet we have verified that off-diagonal instabilities are still present. 
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Xq — 


u.uo 


Fitted moments 




2+6+A 


U.izo ± U.UiU 




n 1 /in -1- n nns 

U.i^tU IC U.UUo 


O + 0+ i 


n 1 Qs -1- n nnn 

u.ioo ± u.uuy 


/+4+0+O 


u.i4z ± U.uuy 


O+0+ r +y 


U.iz4 ± U.UU / 


9 -1-4 -1-^1 -J- 7 
z-t-4-t-d-t- ( 


n 1 AT -1- n nriQ 


3+4+5+6+7 


0.1256 ± 0.0049 


3+4+5+6+8 


0.1247 ± 0.0050 


2+4+5+6+8 


0.1242 ± 0.0042 


2+4+5+7+8 


0.1254 ± 0.0044 



Table 5: 

Fits of as{Mz) from the evolution of an increasing number of moments, with optimal Xq. 

which this happens for more than one experimental point. 

Clearly, the maximum number of moments which may be included in the fit before an off- 
diagonal instability appears will be larger when correlations are lower. This means that first, 
it is not advantageous to further increase the number of scales beyond three (compare with 
Table second, it is convenient not to fit simultaneously neighbouring moments (compare 
with Table |l|); finally, it is convenient to choose a value of Xq which is as low as possible since, 
as discussed above, correlations rapidly increase with xq. Indeed, we find that, with three scales 
and Xq = 0.03, fits of up to five moments are possible, while with xq = 0.1 fits of at most four 
moments are possible without incurring in instabilities. With xq = 0.01 fits of up to seven 
moments are possible, but this choice of xq requires a considerable amount of extrapolation. 

The results for representative fits with xq = 0.03 are shown in Table ^ as an increasing 
number of moments is fitted. It is clear that both the size of the error and the stability of 
the central value improve as the number of fitted moments increases. Stability of the central 
value of as{Mz) is found with the largest number of fitted moments allowed before the onset of 
off-diagonal instabilities. 

4.4 Variation of the fit parameters 

We now turn to the effect of variations of the truncation point Xq and of the range and number 
of scales. Comparing the results of Table |^ with those obtained with a larger value of Xq shows 
that with higher xq it is not possible to achieve satisfactory stability of the best-fit value of ctg 
before the onset of off-diagonal instabilities. Some representative results obtained with xq = 0.1 
are displayed in Table similar results are obtained for other values xq > 0.03. If instead xq 
is lowered to Xq = 0.01 or below, the error on the best-fit as does not improve further, despite 
the fact that fits with a larger number of moments are possible, and remains in fact somewhat 
larger than the best fit of Table ^. This is consistent with the fact that lowering xq below 0.03 
does not introduce any new experimental information. 
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Xq — 


n 1 


Fitted moments 




2+6+A 


U.izo ± U.UUo 


9-U4-l-fi 
z-t-^-t-u 


n 1 on -1- n m n 


3+5+7 


0.121 ± 0.015 


2+4+6+8 


0.137 ± 0.008 


3+5+7+9 


0.140 ± 0.012 


2+4+5+7 


0.114 ± 0.009 



Table 6: 

Fit of as{Mz) from the evolution of an increasing number of moments with large Xq. 



range (GeV^) 


as 


20-70 


0.1242 ± 0.042 


20-100 


0.1239 ± 0.049 


30-70 


0.1239 ± 0.052 


30-100 


0.1249 ± 0.059 



Table 7: 

Dependence of the value of on the range for the 2 + 4 + 5 + 6 + 8 fit with xq = 0.03 



Coming now to scale choices, the first issue is the dependence of our results on the reference 
scale Qq. If we were fitting all moments, results would be entirely independent of this scale, 
except insofar as different choices of Qq correspond to different choices for the first guess of the 
values of the moments in the minimization routine. However, since only a subset of moments is 
fitted, there might be a residual dependence on Qg, due to the fact that the scale dependence 
of the central values of the moments which are not fitted might not agree completely with the 
predicted scale dependence. We have found that, in fact, the choice of initial scale may affect 
the onset of off-diagonal instabilities: for instance, the fit of moments 2 + 4 + 5 + 6 + 8 turns 
unstable if Ql > 40 GeV^. This is not surprising, in view of the fact that the central values 
of low moments are less accurate at this scale, because of the need to extrapolate at small x. 
Indeed, the instability does not appear if evolution is performed using the method of Ref . , 
instead of that of Ref. []TT| : that method is less accurate, but it does not require knowledge of 



low moments since the evolution matrix is triangular (as discussed in section |]) . Nevertheless, 
we find that, in all cases in which a stable fit is obtained, the results turn out to be essentially 
independent of the choice of the reference scale Qq. 

Next, we consider the dependence on the choice of fitting scales. The range of scales cannot 
be widened much without requiring considerable extrapolation. The effect of small variations 
is displayed in Table 0. It is apparent that no significant dependence is found, and in fact the 
smallest error is obtained when 20 < < 70 GeV^. If instead the number of scales is varied, 
we find that with only two scales the quality of the fit deteriorates considerably: for instance, 
the error on as from the fit of moments 2 + 4 + 5 + 6 + 8 with xq = 0.03 goes up to cr = 0.0077 
if only two scales are used. If four or more scales are used, fits with five moments become 
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Figure 3: 

as a function of as(M^) for the fit of moments 2 + 4 + 5 + 6 + 8 with xq = 0.03. 

unstable because of excessive correlations, and it becomes impossible to find a stability region. 
Hence the choice of three scales in the range 20 < < 70 GeV^ appears to be optimal. 

4.5 Best fit 

Having tested how the quality of the fit varies with all the choices enumerated in the previous 
subsections, we select the fit architecture that maximizes the stability and minimizes the error. 
Our conclusion is that a reliable and stable determination of agiMz) is obtained with xq = 0.03, 
five moments (Table |^) and three scales 20 < < 70 GeV^. Specifically, the smallest error is 
obtained with the 'symmetric' combination of moments 2 + 4 + 5 + 6 + 8, which we take as our 
baseline result. Note that the central value of as coincides with that which is obtained from 
single-moment fits (Table |]) when the error is stationary, i.e. for 5 < n < 6. In order to obtain 
a more accurate determination of the statistical error, we have studied the dependence of the 
on the value of ag, displayed in Fig. |^. The is asymmetric about the minimum, rising 
more slowly as ag decreases. We arrive thus at the determination 

as[,Mz) = 0.124 t (stat.) . (4.32) 



19 



4.6 Theoretical uncertainties 



The main sources of theoretical uncertainty are higher-order perturbative and higher-twist cor- 
rections, as well as heavy quark threshold effects. First of all, our fits are based on leading-twist 
evolution of structure functions, so one should worry about possible power corrections. The 



largest corrections of this kind are target mass corrections, which are known analytically p7 



We have checked that these corrections are less than 1% on any of the moments included in 



our fits. Because higher-twist corrections to the operator product expansion are known |^ to 
be significantly smaller than target-mass corrections, we conclude that all power corrections 
are entirely negligible in our analysis, thanks to the relatively high cut > 20 GeV^ which 
we have imposed. Also, one might worry that very high moments might be sensitive to elastic 
contributions to the cross section [OS. We have verified that such contributions are below 1% 



for all moments if > 30 GeV^, and reach at most 3% for the 8-th moment at = 20 GeV^, 
which is the highest fitted moment. We conclude that these contributions are also negligible. 
Nuclear corrections to the deuterium structure function affect the initial values of the fitted 
moments but not their scale dependence, and are thus immaterial, up to nuclear higher twist 
corrections, which are not expected to be significantly larger than standard higher twist terms 
(except possibly at very small x) [pO| . 

We are then left with uncertainties related to higher-order perturbative corrections and to 
heavy quark thresholds. The position of thresholds = kthM^ has been varied in the range 
0.3 < kth < 4. Higher-order corrections have been estimated by varying the renormalization 
scale /i^g„ = krenQ"^ in the standard way |^, with 0.3 < kren < 4. Notice that, because the 



structure function is evolved directly in the DIS scheme |]13|, there is no factorization scale 
dependence. The of fit is almost independent of the choice of kth, while it tends to 
increase considerably if k^en < 0.5 or kren > 2; in particular, with kren < 0.25 we were unable 
to obtain a stable fit. 

The dependence of the value of as(M^) on kth and kren for the fit of moments 2-|-4-|-5-|-6-|-8 
with xq = 0.03 and three scales is displayed in Figure ^. The associated uncertainties are 
estimated to be 

a(thresh.) = t ^^.^ ; a(ren.) = 1 '^.Zl . (4.33) 

The dependence on the position of thresholds is, predictably, very weak, given that the b 
threshold is close to the edge of our range, and falls outside it as soon as kth < 0.8. The 
dependence on renormalization scale turns out to be also reasonably weak. 

4.7 Result and comments 

Our final determination of the strong coupling is 

as{Mz) = 0.124 t (exp.) t S (th-) = 0.124 t (total) . (4.34) 

The error on the result is dominated by statistical uncertainties, consistent with our expectation 
that the method of analysis used here minimizes theoretical uncertainties. 

The value of as has been previously extracted from QCD fits to the BCDMS data, with the 
result [gg] asiMz) = 0.113±0.003(exp.) ±0.004(th.), and to the NMC data, with the result |2B 



as{Mz) = 0.117t 0016 • There are two main differences between these previous determinations 
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Figure 4: 

Dependence of the best-fit value of as{Mz) on the position of heavy quark thresholds (left) 
and renormalization scale (right). The band indicates the overall uncertainty. 

and the present one, eq. ([4.34|) . First, the determinations in Refs. [0,^ were based on global 



QCD fits, and thus required the construction of a parton parametrization, which, as discussed 
in the introduction, might be the source of systematic error and bias. Second, they did not 
include a full treatment of correlated systematics: statistical and systematic errors were added 
in quadrature. A reanalysis of the BCDMS data which did include a treatment of correlated 
systematics found instead as{Mz) = 0.118 ± 0.002 (exp.). Our value appears thus in good 



agreement with other determinations from the same data. While a direct comparison of the 
uncertainties is difficult, we find it interesting that, while the overall uncertainty of our value 
is very close to that of the BCDMS determination |^2[, our analysis gives excellent control on 



theoretical uncertainties: the dominant error is experimental and could be improved upon by 
future experiments. 

It is interesting to observe that the central value we find, though in agreement within error 
with current global averages, is on the high side. As it appears from Table ^, the central value 
of as tends to be higher for higher moments. This is very suggestive of soft gluon resummation 
effects: as is well known |^2|, leading higher order corrections at large are resummed by the 
replacement —>■ Q'^/N in the argument of as in the evolution equation for the A^-th (full) 
moment. This means that the effective value of as is larger for the evolution of higher moments. 
Our results provide some indication of this effect, and suggest that a better determination of as 
could be obtained from these data if this resummation were included, by generalizing the soft 
gluon resummation formalism to the case of truncated moments. In particular, the inclusion of 
soft gluon effects to all orders is likely to reduce the theoretical uncertainty expressed by the 
dependence on the renormalization scale. 
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5 Concluding remarks 



We have presented a determination of from scaling violations aimed at minimizing the 
sources of theoretical bias which might be cause of concern in existing determinations of 
from deep inelastic scattering experiments. This has been accomplished by avoiding the use of 
a parton parametrization: thanks to the use of truncated moments, we have directly fitted the 
scaling violations of a physical observable. Truncated moments in turn have been determined 
by means of a bias-free parametrization of the structure function, inferred from the data, which 
retains all the information on experimental errors and correlations. 

It is interesting to compare our final result for to other determinations from deep inelastic 
scattering, and in particular those obtained using the same data set [^^ -p^. Whereas in our 



determination theoretical errors are small and fully under control, the price to pay for this 
is that the experimental error is larger than in other determinations. It is unclear to which 
extent this is a trade-off, or a consequence of the need to impose restrictive cuts in the Q"^ 
range in order to deal with truncated moments. However, it does suggest that a substantially 
more precise determination could be obtained using data with either smaller statistical error 
(especially for deuteron data), or spanning a wider Q"^ range, or both, such as could be obtained 
for instance at the planned EIC facility ||33|| . 

As far as the central value of as is concerned, it is suggestive that our determination is 
on the high side, since this is what one would expect in the presence of sizable soft gluon 
resummation effects. These in turn are expected to be more important in our determination, 
since the kinematic cuts which we imposed give more weight to the large- a: region, where such 
effects are larger. Therefore, we expect that a resummed version of our analysis might lead to 
a more accurate result for a^, and provide explicit evidence for soft gluon resummation. 

Finally, our analysis provides an explicit demonstration of the power of methods of analysis 
based on the direct determination of the probability density of physical observables from the 
data, and their use coupled with the bias-free computational method which is afforded by the 
use of truncated moments. 
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